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FOREWORD 


The  work  described  in  this  report  was  sponsored  by  the  Independent  Research 
(IR)  Board  of  the  Naval  Surface  Weapons  Center  (NSWC).  It  was  performed  during 
the  period  October  1976  to  September  1979. 

The  traditional  uses  of  radar  have  been  to  locate  a  target  and  to  track 
its  direction  and  velocity.  In  recent  years,  interest  in  a  third  use  of  radar 
has  evolved:  to  identify  an  object  from  the  scattered  return.  The  ultimate 
purpose  of  this  work,  which  is  presently  continuing,  is  to  devise  a  method  for 
identifying  a  complex  object  from  the  return  of  a  radiated  impulse  or  other 
suitably  shaped  short  pulse.  In  order  to  be  able  to  do  this,  the  impulse  re¬ 
sponse  radiation  of  simple  objects  must  be  studied  to  ascertain  the  relation¬ 
ships  between  the  shape  of  an  object,  the  material  of  which  it  is  made,  its 
aspect  angle,  etc.,  and  the  scattered  return. 

This  report  describes  the  NSWC  time-domain  range  returns  from  simple  ob¬ 
jects  and  first  attempts  to  process  them.  Prony’s  algorithm  is  described  and 
efforts  to  improve  it  are  detailed.  The  application  of  Prony's  algorithm  to 
the  returns  is  described  and  the  results  are  shown.  Problems  with  clutter  and 
noise  are  discussed,  and  directions  for  future  efforts  are  recommended. 

The  ultimate  outcome  could  result  in  a  radar  target  identification  system, 
which  would  be  of  use  to  all  the  military  services. 

The  authors  wish  to  gratefully  acknowledge  the  helpful  critiques  and  en¬ 
couragement  of  Dr.  Ronald  J.  Gripshover,  in  particular,  and  of  the  NSWC  IR 
Board,  in  general. 

This  report  has  been  reviewed  by  W.  S.  Orsulak,  Head,  Special  Applications 
Branch,  and  K.  C.  Baile,  Head,  Advanced  Projects  Division. 
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INTRODUCTION 


The  objective  of  this  project  is  to  devise  a  method  for  the  identification 
of  complex  objects,  based  on  the  time-domain  returns  from  a  radiated  short 
(subnanosecond)  pulse.  This  report  details  the  present  status  of  the  project; 
describes  the  time-domain  range  and  experimental  apparatus;  describes  the  pres¬ 
ent  status  of  the  apparatus  and  the  needs  for  refinement  and  improvement.  The 
signal  processing  is  discussed,  and  in  particular,  the  problems  with  Prony’s 
algorithm  are  stated  and  the  improvements  to  the  algorithm  for  coping  with  the 
problems  are  given.  The  results  of  applying  these  to  synthetic  waveforms  are 
given  and  the  areas  of  need  for  further  effort  pointed  out.  An  attempt  to 
apply  the  algorithm  to  experimental  data  is  described,  and  the  results  thereof 
are  discussed.  Finally,  directions  for  further  effort  are  suggested. 


BACKGROUND 


Baum1  expanded  the  solution  of  a  free-radiating  shape  (no  forcing  func¬ 
tion)  in  terms  of  its  singularities,  or  complex  resonances,  in  the  Laplace 
(complex  frequency  or  "s")  domain.  His  expansion,  for  finite-size  objects, 
takes  the  form 


F  (s) 


N  A. 

E  _ 1 

t-i  8  -  si 


(D 


where  the  are  the  complex  resonances  or  poles  of  the  waveform,  and  the 
A^  are  the  coefficients  or  residues.  In  theory,  the  number  of  resonances,  N, 
is  infinity,  but,  in  practice,  the  high-frequency  resonances  have  negligibly 
small  residues  that  render  them  undetectable.  Hence,  a  finite  number  of  reso¬ 
nances  are  considered.  Baum  showed  that  the  poles,  s^,  depend  upon  the  ob¬ 
ject:  its  size,  its  shape,  and  the  material  of  which  it  is  made.  He  also 
showed  that  the  residues,  depend  upon  the  initial  waveform:  its  shape, 
polarization,  and  coupling  to  the  object.  The  A^  also  depend  upon  the  aspect 
angle  of  the  object  in  relation  to  the  transmitter.  Note  that  the  poles,  s^, 
depend  only  upon  the  object  and  are  independent  of  the  transmitted  waveform. 
This  theoretically  means  that  we  can  observe  the  same  set  of  resonances,  inde¬ 
pendent  of  the  direction  the  observer  is  from  the  object. 


If  we  take  the  inverse  Laplace  transform  of  Equation  (1),  we  have 


f  (t) 


N 

E 

i»1 


Hi  esifc 


(2) 


This  is  just  the  sum  of  a  set  of  complex  exponential  terms.  The  poles,  S£, 
have  negative  real  parts.  That  is,  the  free-radiating  waveform  must  damp  out 
in  the  absence  of  new  energy  supplied  to  the  object.  The  function,  f(t),  is 
real,  since  it  is  a  physical  observable. 
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Prony’s  algorithm2  is  a  method  for  extracting  the  poles  and  residues  from 
a  waveform  of  the  form  of  Equation  (2)*  The  aim  is  to  apply  this  algorithm  to 
the  time-domain  radiation  from  an  object  to  extract  its  poles. 


EXPERIMENTATION 


The  data  gathering  operation  is  performed  on  a  9.1-m  (30-ft)  diameter 
ground  plane  shown  in  Figure  1.  The  ground  plane  is  constructed  of  aluminum 
plates  welded  together  to  form  a  continuous  conducting  surface. 


The  electronic  apparatus,  located  beneath  the  ground  plane,  is  shown  in 
Figure  2.  The  generator  is  an  IKOR  R-100  impulse  generator,  which  puts  out  a 
subnanosecond  Gaussian  pulse  with  an  amplitude  of  approximately  1  kV.  The 
repetition  rate  is  approximately  100  Hz.  At  the  output  of  the  impulse  genera¬ 
tor  is  an  FXR  feedthrough  capacitive  pickoff  probe,  the  output  of  which  deliv¬ 
ers  a  trigger  pulse,  through  appropriate  attenuation,  to  a  sampling  oscillo¬ 
scope.  The  signal  must  go  through  a  delay  line  in  order  to  allow  the  sampling 
scope  time  to  initiate  the  trace  on  receipt  of  the  trigger  pulse.  This  delay 
line  consists  of  about  15.2  m  (50  ft)  of  RG-331  cable.  After  the  delay  cable, 
there  is  a  measuring  probe,  which  was  designed  and  built  by  IKOR.  It  is  a 
capacitive  feedthrough  probe  with  an  attenuation  of  49  ±1  dB  and  a  flat  fre¬ 
quency  response  from  5  MHz  to  3  GHz.  After  this  probe,  an  RG-9  cable,  approxi¬ 
mately  O.Q-m  (2.5-ft)  long  carries  the  signal  to  the  base  of  the  transmitter 
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antenna*  With  the  exception  of  the  output  of  the  impulse  generator  and  the  FXR 
pickoff  probe/  which  uses  Type  N  connectors,  the  connectors  from  the  generator 
to  the  base  of  the  transmitter  antenna  are  all  of  type  HN.  This  type  of  con¬ 
nector  was  chosen  in  anticipation  of  future  higher  voltages. 


The  transmitter  antenna  is  a  5.2-m  (17-ft)  wire  monopole.  Its  upper  end 
is  supported  from  the  ceiling.  The  pulse  is  radiated  by  the  monopole  and  re¬ 
ceived  by  the  receiving  antenna.  Meanwhile,  the  radiated  pulse  strikes  the 
object,  giving  a  return  from  the  object  that  arrives  at  the  receiving  antenna 
a  few  nanoseconds  after  the  initial  pulse. 

The  receiving  antenna  is  a  half  TEM  horn,  0.3-m  (1-ft)  long  connected  to 
the  stub  of  a  Type  N  wall  connector  mounted  in  the  ground  plane.  The  other 
half  of  the  TEM  horn  is  imaged  by  the  ground  plane.  This  simulates  a  TEM  horn 
with  an  angle  of  16*.  These  particular  antennas  (monopole  transmitter  and  TEM 
horn  receiver)  were  chosen  as  a  start  because  the  monopole  transmits  with  high 
fidelity  the  source  waveform,  while  the  TEM  horn  receives  and  reproduces  with 
high  fidelity  the  waveform  of  the  electric  field  of  the  radiation.  Reference  3 
gives  the  details  of  various  antenna  configurations  as  transmitters  and 
receivers. 
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Figure  3  shows  the  relative  positioning  of  the  transmitter  (monopole) , 
the  receiver  (TEM  horn) ,  and  the  object.  (The  paint  can  shows  the  position 
of  the  object.)  They  were  positioned  in  line  so  as  to  simulate  a  monostatic 
system;  i.e.,  the  object  is  in  the  same  direction  from  both  transmitter  and 
receiver.  A  check  revealed  that  the  presence  of  the  transmitter  between  the 
receiver  and  object  does  not  affect  returns  within  the  sensitivity  of  the  in¬ 
strumentation.  The  distance  between  the  monopole  and  the  TEM  horn  is  2.3  m 
(7.5  ft).  The  reason  for  this  distance  is  to  protect  the  receiving  equipment 
from  the  amplitude  of  the  transmitted  pulse,  until  appropriate  duplexing  equip¬ 
ment  can  be  built.  The  part  of  the  object  nearest  the  transmitter  is  0.6  m  (2  ft) 
from  the  transmitter.  This  will  assure  a  delay  of  about  4  nsec  from  the  time 
the  radiated  pulse  is  received  by  the  TEM  horn  until  the  return  from  the  object 
arrives.  (Light  travels  30  cm  in  1  nsec) •  This  4-nsec  delay  is  necessary  to 
allow  for  the  trail-off  from  the  transmitted  pulse. 


Figure  3.  Relative  Positioning  of  Transmitter  (center).  Receiver 
(upper  left) ,  and  Object  (lower  right) 


Upon  receipt  by  the  TEM  horn,  the  signal  is  fed  to  a  sampling  oscilloscope 
through  an  RG-9  cable  and  appropriate  attenuators.  The  sampling  oscilloscope 
is  a  Hewlett-Packard  Model  181A  dual-beam  oscilloscope.  It  has  a  Model  181 1A 
sampling  time  base  and  vertical  amplifier  plug-in.  The  sampling  head  is  a 
Model  1430C.  All  connectors  in  the  receive  circuit  are  of  Type  N. 
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On  the  back  of  the  sampling  scope's  main  frame  are  outputs  for  monitoring 
waveforms  by  external  recording  equipment.  The  waveforms  that  can  be  monitored 
are  the  two  traces  displayed  on  the  screen  and  the  horizontal  sweep  (sawtooth) 
waveform.  A  Hewlett-Packard  3960  analogue  tape  recorder  was  used  to  record  the 
waveforms  for  later  data  reduction  on  a  PDP  1145  computer.  The  idea  was  to 
record  about  30  traces  of  the  returns  from  each  object  and  later  to  have  the 
computer  digitize  the  returns  and  compute  an  average  for  each  return.  A  major 
advantage  of  using  an  analogue  tape  recorder  is  that  the  original  traces  would 
be  recorded  permanently  and  could  be  referred  to  as  often  as  needed.  However, 
the  tape  recorder,  even  within  the  manufacturer's  specification  and  with  the 
use  of  filters,  puts  out  noise  at  such  a  level  as  to  render  proper  digitizing 
and  averaging  of  data  extremely  difficult  and  uncertain.  Therefore,  to  obtain 
results  for  this  report,  traces  of  returns  were  enlarged  to  30.5  cm  (12  in.)  by 
38.1  cm  (15  in.)  and  hand  digitized.  More  will  be  said  about  this  process 
later. 

The  objects  whose  returns  are  being  studied  are  a  hemisphere,  a  cone,  and 
a  cylinder,  which  are  shown  in  Figure  4.  They  are  each  12.7  cm  (5  in.)  in 
diameter  at  the  base.  The  cone  has  a  30°  angle,  and  the  cylinder  is  25.4  cm 
(10  in.)  in  height.  They  are  all  solid  aluminum. 
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Figure  5  shows  sampling  oscilloscope  traces  of  the  impulse  generator  pulse* 
Figure  5a  shows  it  as  monitored  by  the  IKOR  probe  at  the  output  of  the  delay 
line  and  just  before  transmission  by  the  monopole  antenna.  Figure  5b  shows  the 
pulse  as  received  by  the  TEM  horn  receiving  antenna.  The  wiggles  in  the  trail- 
offs  of  these  traces  arise  from  reflections  caused  by  connectors  in  the  trans¬ 
mission  lines  to  and  from  the  antennas. 

Figure  6  shows  sampling  oscilloscope  traces  of  waveforms  in  the  transmit¬ 
ter  (IKOR  probe)  and  receiver  (TEM  horn)  circuits,  with  no  object  on  the  ground 
plane.  The  lower  trace  shows  the  waveforms  in  the  transmitter  circuit.  The 
first  pulse  in  the  lower  trace  is  the  initial  impulse  generator  pulse  as  seen 
in  Figure  5a.  The  second  pulse  is  the  reflection  of  the  initial  pulse  at  the 
base  of  the  transmitter  monopole.  Notice  the  third  pulse  that  appears  about 
40  nsec  after  the  initial  pulse.  It  represents  a  reflection  of  energy  from  the 
tip  of  the  transmitter  monopole,  a  distance  of  r.2  m  (17  ft)  from  the  base. 

The  upper  trace  shows  the  waveforms  in  the  receiver  circuit.  The  first  pulse 
is  the  impulse  generator  pulse.  After  its  trailoff,  except  for  clutter,  there 
is  a  clear  time  window  of  about  20  nsec.  The  first  small  negative  peak,  just 
after  the  second  vertical  grid  line  to  the  right  from  the  center,  is  a  reflec¬ 
tion  from  the  nearest  wall  of  the  building.  The  large  oscillations  immediately 
after  that  are  caused  by  ringing  by  the  transmitter  monopole  antenna  and  re¬ 
flections  from  the  ground  plane  edges  and  building  walls.  Figure  7  shows  the 
same  waveforms  as  Figure  6  except  that  the  aluminum  cylinder  is  placed  in  the 
object  position.  No  changes  are  seen  in  the  lower  trace  (the  transmitter  cir¬ 
cuit  waveform) .  The  upper  trace  (the  receiver  circuit  waveform)  is  the  same  as 
before,  except  that  the  return  from  the  cylinder  is  seen,  immediately  after  the 
impulse  generator  pulse  trailoff.  The  pictures  of  traces  from  target  returns 
will  show  the  10-nsec-time  window  from  the  third  (from  the  left)  to  the  fifth 
(center)  grid  line  of  Figure  7. 

Figure  8  shows  expanded  traces  of  returns  from  the  cylinder,  the  cone,  and 
the  hemisphere.  The  trace  marked  '’clutter*  is  what  is  seen  in  the  same  time- 
frame  when  no  object  is  present.  These  photographs  were  enlarged  to  30.5  cm 
(12  in.)  by  38.1  cm  (15  in.)  in  such  a  way  as  to  make  the  squares  of  the  grids 
the  exact  same  size.  Xerox  copies  of  the  enlargements  were  made  on  thin  paper  • 
These  copies  were  carefully  checked  for  distortion  by  measuring  the  squares  on 
the  grid  and  comparing  with  the  squares  on  the  photographed  enlargements.  The 
worst  distortion  was  found  to  occur  near  the  edges,  away  from  the  important 
parts  of  the  traces,  which  showed  negligible  distortion.  The  curves  were  then 
hand  digitized  from  the  thin-paper  copies  on  an  Elographics  digitizing  table, 
Every  effort  was  made  for  even  spacing  along  the  time  axis.  Approximately  201 
points  were  digitized  from  each  curve.  A  cubic  spline  routine  similar  to  that 
described  in  Reference  6  was  used  to  smooth  each  data  set  and  interpolate  to 
512  points,  equally  spaced  along  the  time  axis.  The  smoothed  data  sets  were 
submitted  to  our  CDC  6700  computer  for  processing. 
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BEFORE  TRANSMISSION 


200  mv/div 
1  ns/div 


RECEIVED  BY  TEM  HORN 


Figure  5.  Impulse  Generator  Pulse 


Figure  6.  Transmitted  and  Received  Waveforms;  No  Object  on  Ground  Plane 
(5  nsec/div;  200  mV/div)  (Lower  trace — transmitted  pulse  and  its 
reflections*  Upper  trace — as  received  by  TEM  Horn) 


Figure  7.  Transmitted  and  Received  Waveforms;  Cylinder  on  Ground  Plane 
in  Object  Position  (5  nsec/div;  200  mV/div)  (Lower  trace — transmitted 
pulse  and  its  reflections.  Upper  trace — as  received  by  the  TEM  horn) 


PRONY'S  ALGORITHM 


This  section  gives  the  exact  method  used  by  Prony's  algorithm  for  deter¬ 
mining  the  poles  and  residues  of  a  waveform  described  by  Equation  (2) •  The 
procedure  is  basically  that  outlined  in  Reference  5. 

Let  us  sample  f(t)  at  N  points,  equally  spaced  in  time.  Let  us  set  the 
sampling  times  as  follows: 

t^  *  (k— 1)6,  k  *  1*  2,  •••$  N 


where  6  is  the  time  interval  between  samples.  Also,  let  us  rewrite  Equation  (2) 
as 


s  f(tk)  -  ^  V 


a»(k-')8  "  .  „  k-i 


m*1 


m  m 


(3) 


where 


xm  »  e 
m 


s  6 
m 


Form  the  following  set  of  equations: 


A1 

+  a2 

+  ...  + 

*  f1 

Vi 

+  a2x2 

+  ...  +  AjjXjj 

“  f2 

(4) 

Vi2 

+  AjXj 

+  +  Vm 

*  f3 

Ml"  1  +  A2x2N  1  +  *•*  +  AMxMN  1  "  fN 
Let  ^  represent  the  roots  of  an  Mth  order  polynomial 


M+1 


Z  cmx 
m=1  m 


m- 1 


( x-x  ^ )  ( x— X2 )  •  • .  ( x— Xjyj)  *  II^Cx-x^ 


(5) 


where 


CM+1  *  1 

Let  us  multiply  the  first  equation  of  (4)  by  ct,  the  second  equation  of  (4)  by 
C2 1  etc.,  for  the  first  M+1  equations.  Then,  using  Equation  (5),  we  fora  the 
equation 
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f1c1  +  f2c2  + 


•••  +  fM+1  “  0 

Repeat  this  operation,  beginning  with  the  second  equation  of  (4),  thereby  forming 
f 2C 1  +  f3C2  +  +  fM+2  =  °* 

Repeat  this  operation  until  M  such  equations  are  formed 
f1c1  +  f2C2  +  "*  +  fM°M  =  -fM+1 

f2c1  +  f3C2  +  •••  +  fM+i  CM  “  “fM+2  (6) 


fMC1  +  fM+1C2  +  •••  +  f 2M-1  CM  =  ~f2M 
We  may  express  Equation  (6)  in  matrix  form  as  follows 
P  C  *  B , 

where 


the  data  sample  matrix, 


the  vector  of  the  coefficients,  and 
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(7) 


Equation  (6)  is  solved  for  the  c*s.  Then,  we  use  the  c's  in  Equation  (5)  to 
solve  for  the  x's,  from  which  we  can  obtain  sm  using 

sm  =  — 1—  In  (8) 

m  6 

The  A's  can  be  found  via  Equation  (4). 

The  foregoing  discussion  detailed  the  exact  procedure  for  calculating 
the  poles  and  residues  using  Prony*s  Algorithm.  If  M  is  the  number  of  poles 
to  be  determined,  we  need  2M  sample  data  points  for  the  exact  calculation. 
Often,  there  are  more  than  2M  data  points  available.  When  more  than  2M  points 
are  used,  the  problem  is  overde ter mined.  The  next  section  will  describe  new 
procedures  for  handling  the  overde ter mined  case. 

In  activating  the  algorithm  as  described,  one  must  know  in  advance  the 
number  of  poles  or  guess,  if  the  number  is  not  known.  If  the  estimated  number 
is  too  low,  the  results  will  be  nonsense.  If  the  estimate  is  too  high,  the 
algorithm  will  give  the  correct  poles,  plus  some  extra  (curve-fitting)  poles 
with  low  residues.  However,  to  obtain  meaningful  results,  the  estimate  must 
not  be  too  far  above  the  correct  number.  ANALYSIS  TECHNIQUES  will  describe 
what  has  been  done  to  handle  these  shortcomings. 


ANALYSIS  TECHNIQUES 


The  most  critical  procedure  in  Prony*s  method  (in  the  overde ter mined  case) 
is  finding  the  least-squares  solution  of 

F  C  =  B  (9) 

for  £,  the  vector  of  coefficients  for  the  polynomial.  The  linear  least- 
squares  procedure  for  solving  Equation  (9)  for  C  is  given  in  Reference  6. 

Once  the  coefficients  are  accurately  known,  finding  the  roots  of  the  polynomial 
is  generally  straightforward.  The  residues  are  less  important  for  most  pur¬ 
poses,  because  they  are  not  necessarily  characteristic  of  the  object  being 
irradiated. 


GROUT  ROUTINE 

The  original  SEMPEX  program  used  the  CROUT7  method  for  finding  the  polyno¬ 
mial  coefficients.  In  the  CROUT  subroutine,  Equation  (9)  is  multiplied  by  FT 
to  obtain 

FT  F  C  =  FTB  (10) 
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(11) 


T 

Here,  F  is  the  transpose  of  F.  The  vector, 

C*  »  (FTP)-1FTB 

which  is  the  unique  solution  of  Equation  (10)  is  the  least-squares  solution 
of  Equation  (9).  The  CROUT  method  uses  a  modification  of  the  Gauss  Reduction7 
to  find  the  solution.  It  does  not  work  well  when  the  columns  of  F  are  close  to 
linear  dependence  because  the  determinant  of  FTF  will  be  small  and  the  prob¬ 
lem  nearly  singular. 


HFTI  ROUTINE 

We  found  the  HFTI8  routine,  available  in  the  Naval  Surface  Weapons  Center 
(NSWC)  Library  of  Mathematics  Routines,  to  be  more  suitable  for  our  purposes. 

It  is  more  accurate  because  most  inner  products  are  computed  in  double  preci¬ 
sion.  It  is  able  to  deal  with  the  problem  of  rank  deficiency  when  the  columns 
of  F  are  close  to  being  linearly  dependent.  The  HFTI  routine  can  also  be  used 
to  get  an  indication  of  the  number  of  poles  present  in  the  signal.  In  the 
process  of  computing  the  least-squares  solution,  the  HFTI  routine  transforms 
the  matrix,  F,  to  upper  triangular  form  with  an  orthogonal  matrix,  Q,  and  an 
orthogonal  permutation  matrix,  P,  so  that  QFP  »  R,  and  R  has  zeroes  below  the 
diagonal.  The  matrix,  P,  reorders  the  columns  so  that  the  absolute  values  of 
the  diagonal  elements  are  monotonically  decreasing.  The  HFTI  routine  then 
calculates  a  condition  number  C|  for  this  triangular ized  data  matrix  by  di¬ 
viding  the  maximum  (upper  left)  diagonal  element,  by  the  minimum  (lower 

right)  diagonal  element,  a^.  That  is,  *  ail/aMM*  The  rank  k  of  F  is  the 

number  of  linearly  independent  columns  of  F.  If  the  number  of  columns  M  is 
greater  than  k,  the  last  M-k  diagonal  elements  of  R  will  be  0.  Ttoat  is,  as  i 
is  increased,  a^  decreases  raonotomically  and  goes  abruptly  to  zero  when 
i>k.  In  this  case,  a  90  for  M>k#  In  fche  presence  of  noise,  there 

is  still  an  abrupt  decrease  in  a^  (and  an  abrupt  increase  in  C|)  as  i  becomes 
greater  than  k,  although  a^  is  no  longer  equal  to  zero  (and  C|<»)  for  i>k.  As 
the  noise  is  increased,  the  decrease  in  a^  (and  increase  in  C|)  as  i  becomes 

greater  than  k,  is  less  abrupt. 

Assuming  no  noise  in  the  data,  the  number  of  linearly  independent  columns 
(or  the  rank  k)  of  F  is  usually  equal  to  the  number  of  poles.  This  is  true  for 
the  following  reason.  Assume  m  is  the  actual  number  of  poles,  and  M  is  the 
assumed  number  of  poles  with  M>m.  Let  N  be  the  number  of  data  points.  Now  F 
is  defined  to  be 


fN-M  fN-M+1  fN-1 
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If  there  is  no  noise ,  then 


f:  * 


E  A»x 
<=1  K  K 


(12) 


where 


l  C_  Xkn_1  +  xk  =  0 
n=1  n  K  K 


Denote  the 


j  'th  column  of  F  by  Fj .  By  Equation  (12) 

r—  — . 


m  . 

£  At  x? 
k=1  K  K 


in  .  .  . 

E  Av  xP+ 1 
k=1  K  K 


m 


L,  Ak  xic 


j+N-M-1 


k=1 


By  Equation  (13) 


£  C  F  +  F  =0. 
n=1  n~n  -m+1  “ 


(13) 


(14) 


Thus,  the  number  of  linearly  independent  columns  cannot  be  greater  than  m. 

It  is  possible  that  there  could  be  fewer  than  m,  but  this  would  be  very 
unusual. 

Thus,  if  the  number  of  linearly  independent  columns  of  F  is  m,  the  rank 
of  F  is  m  and  the  last  M-m  diagonal  elements  of  R  would  be  0  in  the  absence 
of  noise.  We  could  then  conclude  that  the  number  of  poles  is  probably  m  and 
restart  the  algorithm  assuming  m  poles  instead  of  M. 

In  practical  work,  we  check  to  see  when  M  becomes  greater  than  m  by  moni¬ 
toring  the  condition  number  of  matrix  R.  In  successive  runs,  we  note  the 
value  of  M  when  undergoes  an  abrupt  increase.  The  rank  m  is  then  M-1,  the 
value  of  M  just  before  C#  increases  suddenly.  When  is  large  (usually 

greater  than  100  in  practical  situations),  we  say  that  the  data  matrix  is  ill- 
conditioned. 


If  the  signal  is  contaminated  by  noise,  as  is  usually  the  case,  the  situa¬ 
tion  is  somewhat  different.  If  the  noise  is  random,  uncorrelated  noise,  the 
columns  of  F  will  be  linearly  independent  because  the  noise  samples  are  inde¬ 
pendent  of  one  another.  Thus,  the  rank  of  F  will  be  M  and  none  of  the  diagonal 
elements  of  R  will  be  0.  However,  if  the  noise  level  is  relatively  small,  the 
last  M-m  diagonal  elements  will  be  much  smaller  than  the  other  diagonal  ele¬ 
ments.  How  large  then  must  one  set  the  tolerance  for  the  diagonal  elements  to 
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lv 


'C 


determine  how  many  poles  are  actually  present?  Note  that  the  matrix,  Q ,  be¬ 
cause  of  its  orthogonality,  does  not  alter  the  norm  of  the  columns  of  F.  lfce 
matrix,  P,  merely  reorders  the  columns  where  necessary,  so  it  can  be  ignored 
for  this  dicussion. 

Assume  the  diagonal  element  in  column  j  would  be  0,  if  no  noise  were  pres¬ 
ent.  Note  that  Q  will  be  somewhat  different  when  noise  is  present.  Because  of 
the  noise,  the  k'th  element  of  will  be  +  Nj+k-1  ,fhere  is  the 

noise*  If  the  noise  is  random  with  mean  0  and  variance  o^,  then 


N-M  N-M  ^  9 

E(Jc£1  <fj+k-1  +  Nj+k-1>  |  =  kE,  f j+k~1  +  (n-m>°N 


where  E  denotes  expected  value.  The  effect  of  Q  will  be  to  redistribute  this 
quantity  in  the  column  QF\  so  that  all  elements  below  the  j*th  or  diagonal 
element  will  be  0.  Let  us  write  the  kfth  element  of  QFj  as  r^j  +  e^j  where  r^j 

is  the  value  the  element  would  have,  if  there  were  no  noise*  Because  of  con¬ 
servation  of  the  norm  by  Q, 

j  2  N-M 

l  rj:  •  »  I  (16) 

k=1  *3  k=l  1 


with  rjj  =  0,  and 


2  <rkj  +  ek j >  J  “  fj+k + 


These  two  equations  somewhat  justify  the  approximation 


j  .  E  a  (N-M)O^ 


fe(ef.)  =  o (e  .  .  )  = 


Assuming  j  *  M,  then 


0  <ej j  >  - 


If  the  diagonal  element  of  the  last  column  of  QFP  is  smaller  than  this,  the 
number  of  poles  is  probably  less  than  M* 
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NONLINEAR  LEAST  SQUARES 

Contrary  to  what  one  might  assume,  the  normal  linear  least-squares  method 
does  not  give  the  best  fit  to  the  data.  Written  in  terms  of  the  estimated 
solution  for  the  polynomial  coefficients  6  and  the  error  of  the  fit  d,  Prony's 
equation  is 

F  C  =  B  +  d  (21) 

The  usual  linear  least-squares  solution  minimizes  the  norm  of  d.  This  is  the 
error  in  B,  which  does  not  involve  all  of  the  data  values.  Also,  there  are 
actually  2M  parameters  of  the  problem:  the  N  poles  and  the  M  residues.  The 
linear  least-squares  method  optimizes  only  M  parameters  simultaneously. 

The  nonlinear  least-squares  method  of  Evans  and  Fischl4  avoids  these  draw¬ 
backs  at  the  expense  of  a  great  deal  of  extra  computation.  It  is  based  on 
Equation  (21)  and  on 

I*  F1  A  +  e 


with  f 


and  F. 


N-1  N-1 


Here,  is  the  error  vector  involving  the  error  in  fitting  all  the  data  to  be 
minimized  by  a  simultaneous  optimization  of  C  and  A.  Equation  (21)  can  be 
rewritten 

T  A 

B  (C)  f  =  d  (23) 


v 


by  defining  the  N  x  (N-M)  matrix 


There  is  a  corresponding  equation  for  noise-free  data  h^  and  the  correct  values 
of  c  and  no  error# 

BT  <£>  Jl  =  0  ( 24) 

It  follows  then  that 

BT  e  =  BT  (_f-h)  =  d  (25) 

The  problem  is  given  to  find  e^  There  is  no  unique  solution  to  this  as  it 
is  an  underdetermined  problem  because  BT  has  more  columns  than  rows.  This 
represents  a  set  of  N-M  equations  in  N  unknowns.  Another  way  to  look  at  this 
is  that  d  is  a  function  of  the  M  elements  of  C^,  while  e  is  a  function  of  C  and 
additionally  of  the  M  elements  of  A.  The  solution  can  be  made  unique  by 
constraining  A  to  be  a  function  of  C •  The  optimal  way  to  do  this  is,  for  a 
given  _C,  define  A*(C)  to  be  that  which  minimizes  the  norm  of  e^£,  A*)#  Then 
there  will  be  a  matrix  W(£)  such  that 


e(C,  A*)  =  W(C)  d(C)  (26) 

Evans  and  Fischl  show  that  W(C)  is  given  by 

W(C)  =  B(C)  [bT(C)  B(C)J  (27) 
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The  strategy  for  the  nonlinear  least-squares  method  is  to  solve  iterative¬ 
ly  for  the  value  of  C  that  will  minimize  £(£/A*)  in  Bguation  (26).  The  algo¬ 
rithm  begins  by  solving  the  linear  least-squares  problem  of  Equation  (21).  Hie 
solution  of  this  is 

£<°)  =  (ftf)“,1FTB  (28) 

Then,  we  find 

w(0)(^(0))  =  B(C(0))  [  BT(C°)B(C0)]  _1  (29) 

Define  C  * 1  ^  as  the  vector  that  minimizes  the  norm  of  e.  in  the  equation 

e(C)  =  W(0,d(C)  =  W(0)(F  C  -  B)  (30) 

The  solution  will  be 

C*1)  =  [fV0)TW(0)f]-1  FTW(0)T  W(0\b  (31) 

This  solution  for  ^  can  then  be  used  to  compute  ^  ^ ) ,  which  can  be 

used  in  the  least-squares  equation 

e(C)  =  W{ 1  * (F  C  -  B)  (32) 

(2) 

to  find  an  improved  estimate  C'  ' •  This  procedure  can  be  reiterated  as  many 
times  as  accuracy  requires  and  funds  allow.  Evans  and  Fischl  follow  this 
procedure  with  a  gradient  method  to  further  refine  the  solution.  This  was 
introduced  primarily  for  mathematical  completeness  and  is  generally  not  neces¬ 
sary  in  practical  applications.  It  was  not  included  in  our  application  of  the 
algorithm. 


KALMAN  FILTER 

An  alternative  to  the  nonlinear  least-squares  technique  is  a  sequential 
estimation  method  based  on  Kalman  Filter  Hieory.9  It  has  apparently  never 
been  applied  to  this  problem  before  but  shows  considerable  promise  for  situa¬ 
tions  where  multiple  independent  returns  are  available.  In  these  situations, 
it  can  exceed  the  accuracy  of  the  nonlinear  least-squares  method  with  computer 
memory  requirement  comparable  to  linear  least  squares.  It  achieves  this  by 
sequentially  processing  the  data  in  small  blocks,  in  trie  end  obtaining  the 
equivalent  of  one  simultaneous  solution  for  all  the  data  sets.  Hie  data  sets 
can  involve  signals  with  different  residues  and  different  data  spacing  as  long 
as  the  poles  are  the  same. 

The  way  this  method  works  is  as  follows.  Assume  we  have  K  subsets  of 
data,  each  of  which  satisfies  Prony's  equation 

((N,-M)xM)  (M)  (N.-M)  (N^-M)  (dimensions  of  matrices  and  vectors 

3  33  below) 

F(^  C  *  B(3)  +  d(3)  j  -  1,  2,...,  K  (33) 
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Here,  F^)  is  the  j*th  data  matrix  of  dimension  N^-M  by  M,  C  is  the  vector  of 

the  polynomial  coefficients  of  dimension  M,  is  the*^1  th-S  vector  of 

dimension  N^-M  and  is  the  error  vector  also  of  dimension  N^-M.  Now  assume 


that  all  the  data  contain  the  same  set  of  poles,  and,  thus,  the  polynomial 
coefficients  C  should  ideally  be  the  same.  The  object  is  to  obtain,  using  all 
the, data  sets,  the  solution  for  C  that  will  give  the  smallest  value  of 


ill*  ,Jf’ 


i.e.,  the  overall  unweighted  (linear)  least-square  error. 


It  should  be  noted  here  that  the  theory  allows  for  the  minimization  of  a 
weighted  sum  of  squared  errors.  To  obtain  the  minimum  variance  solution,  what 

must  be  minimized  is  t)TW  D  where  D  is  a  column  vector  made  up  of  all  the  (i^'s 

and  W"1  »  E(D^D) ,  the  expectation  value  of  DTD.  In  our  application,  however, 

the  errors  are  assumed  to  be  independent  and  to  have  equal  variance.  Therefore, 
W  becomes  a  constant  times  the  identity  matrix  and  has  no  effect  on  the  theory. 

(V\ 

The  k*th  estimate  for  C,  denoted  Cx  represents  the  simultaneous  least- 
squares  solution  of  the  first  k  subsets  of  data.  TO  continue  the  process  and 
find  ^  also  requires  knowing  which  is  the  covariance  matrix  of  the 

current  estimate  C^.  According  to  the  Kalman  theory,  the  solution  is 


c(k+1)  „  c(k)  +  p(k)p(k+l)T(I+p(k+1)p(k)p(k+1)T)”1(B(k+l)  _  p(k+l)c(k), 
p(k+l )  „  p(k)  _  p(k)F(k+l)T(I+p(k+1)p(k)p(k+1)T)-l  p(k+l)p(k) 


Note  that  the  inverse  is  of  dimension  (N^-MJXCN^-M) . 


The  iteration  begins  by  computing 

pd)  .  (P<»)Trn))-i  <36> 

C{1)  -  P(1)F1TB(1)  (37) 

In  theory  the  number  of  data  points  could  be  different  for  each  data  sub¬ 
set  and  could  be  any  number  >2M.  In  our  application,  we  use  N^  *  2M  points 

for  every  data  set.  This  simplifies  the  computation  somewhat  and  reduces  the 
dimension  of  the  inverse  to  M  x  M.  In  comparison,  the  nonlinear  least-squares 
method  involves  inverses  of  dimension  (N-M)  x  (N-M).  Thus,  for  example,  if 
M  *  14  and  N  *  84,  the  dimension  of  the  nonlinear  least  squares  inverses  would 
be  70  x  70,  involving  4900  elements,  whereas  the  Kalman  method  would  do  this  in 
three  steps  with  each  step  involving  inverses  of  dimension  14  x  14,  or  196 
elements  per  step. 
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APPLICATIONS 


APPLICATION  TO  SYNTHETIC  DATA 


All  of  the  above  methods  work  well,  if  the  noise  level  is  low  and  if  cer¬ 
tain  criteria  on  data  spacing  are  met.  We  have  found  by  experience  that  the 

data  spacing  At  must  be  small  enough  to  satisfy  the  Nyquist  criterion  that 

be  greater  than  or  equal  to  the  highest  frequency  present.  Another  finding  is 
that  the  total  data  span,  T,  must  be  large  enough  so  that  is  less  than  or 

equal  to  the  lowest  frequency  present.  These  criteria  both  seem  reasonable  on 
theoretical  grounds.  There  have  been  cases  where  the  method  worked  without 
satisfying  one  of  these  criteria  but  usually  it  does  not.  A  large  condition 
number  is  often  an  indication  of  trouble  in  this  area. 

The  true  test  of  these  methods  is  their  performance  in  the  presence  of 
noise.  In  order  to  compare  them,  we  have  used  a  set  of  six  pairs  of  poles 
published  by  Van  Blaricum  and  Mittra.10  In  this  report,  they  present  a  method 
that  works  on  synthetic  data  generated  from  these  poles  using  residue  values  of 
1.0  in  the  presence  of  additive  random  Gaussian  noise  with  a  standard  deviation 
0  as  large  as  =  0.01.  With  our  methods,  we  are  able  to  get  reasonable 
results  with  =  0.1. 


The  criterion  we  used  to  measure  the  effectiveness  of  the  above  described 
methods  in  noise  is  the  root  of  the  sum  of  the  squares  (RSS)  of  the  differences 
between  our  calculated  and  Van  Blaricum's  poles.  That  is, 


where  are  Van  Blaricum* s  poles  and  s^  are  the  poles  as  determined  by  our 
methods.  The  lower  RSS  is,  the  more  effective  the  method. 

Table  1  shows  results  obtained  with  the  three  methods  on  synthetic  data 
generated  from  Van  Blaricum* s  poles,  using  the  residues  shown  in  Table  2. 

Random  uncorrelated  Gaussian  noise  was  added  with  a  of  0.1.  Each  synthetic 

return  had  an  independent  noise  sample;  thus,  cases  1  and  2,  which  had  the  same 
poles  and  residues  gave  different  results.  Data  spacing  At  was  0.25  nsec. 

The  Kalman  method  was  used  to  analyze  the  data  in  the  following  manner: 

Iteration  Data  Sets  Used 


1  1 

2  1  and  2 

3  1r  2,  and  3 

4  1,2,3,  and  4 

5  1,  2,  3,  4,  and  5 
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Results  improved  with  each  iteration,  as  more  data  sets  were  included.  In 
these  calculations,  the  Kalman  method  used  150  points  from  each  return  (data 
set)  while  the  other  methods  used  only  80  points  each. 


Table  1.  Results  from  Synthetic  Data  (Van  Blaricum  and  Mittra) 


MS  *V£ih'  4i|  2 


Case  No. /Method 

Linear-Least 

Squares 

Non  li  near -Least 
Squares 

Kalman 

1 

0.0302 

0.0261 

0.0543 

2 

0.0382 

0.0279 

0.0438 

3 

0.0116 

0.0097 

0.0262 

4 

0.0246 

0.0234 

0.0155 

5 

0.0385 

0.0391 

0.0106 

Table  2.  Poles  and  Residues 


_ Residues _ 

Poles/CaBe  No.  _ 1 _  _ 2 _  _ 3  _ 4 _  _ 5 


-0.082 

± 

j 

0.926 

-0.70 

+ 

j 

0.70 

-0.70 

± 

j 

0.70 

0.70 

± 

j 

1.0 

0.70 

+ 

3 

0.70 

-1.0  : 

*  j 

0 

1-0 

-0,147 

± 

j 

2.874 

-0.70 

t 

j 

0.70 

-0.70 

+ 

j 

0.70 

0.70 

± 

j 

1.0 

-0.70 

± 

3 

0.70 

0.70 

j 

0.70 

-0.188 

± 

j 

4.835 

-0.70 

± 

j 

0.70 

-0.70 

± 

j 

0.70 

0.70 

± 

j 

1.0 

0.70 

3 

0.70 

0.70 

+ 

j 

0.70 

-0.220 

± 

j 

6.800 

0.70 

+ 

j 

0.70 

0.70 

♦ 

j 

0.70 

0.00 

± 

i 

1.0 

-0.70 

± 

3 

0.70 

1.0 

± 

j 

0.0 

-0.247 

t 

j 

8.767 

0.70 

± 

i 

0.70 

0.70 

t 

j 

0.70 

0.00 

± 

j 

1.0 

0.70 

♦ 

3 

0.70 

0.0 

j 

1.0 

-0.270 

t 

j 

10.733 

0.70 

± 

i 

0.70 

0.70 

± 

j 

0.70 

0.00 

t 

j 

1.0 

-0.70 

± 

3 

0.70 

-0.70 

± 

j 

0.70 

It  can  be  seen  that  all  three  methods  performed  rather  well.  The  non¬ 
linear  least  squares  generally  outperformed  the  linear  least  squares.  The 
Kalman  method  outperformed  everything  else  except  the  result  of  the  nonlinear 
least  squares  method  on  data  set  3.  To  some  extent  though,  the  results  of 
the  Kalman  method  were  disappointing.  Using  just  the  Case  1  data,  it  should 


have  equaled  the  linear  least  squares ,  but  it  did  not.  Using  just  84  data 
points  in  the  Kalman  method  improved  Case  1  results  to  an  RSS  error  of  0.0407 
but  degraded  the  other  results.  We  do  not  understand  why  this  is  the  case, 
but  given  more  time,  it  would  be  interesting  to  investigate  this  further.  The 
Kalman  method  seems  to  be  programmed  correctly.  We  did  the  runs  in  reverse 
sequence,  cases  5,  4,  3,  2,  and  1  and  got  the  same  final  results.  Also,  for  a 
low-noise  case,  oN  *  0.005,  the  Kalman  and  linear  least-squares  methods  gave 

the  same  results.  Perhaps  the  Kalman  and  linear  least-squares  methods  are  not 
mathematically  equivalent  in  the  presence  of  noise. 

We  found  a  curious  result  involving  the  data  spacing  with  these  computa-  , 
tions.  The  poles  were  first  extracted  using  a  data  spacing  of  At  =  0.2857  nsec 
giving  a  Nyquist  frequency  of  1.75  GHz,  which  exceeds  the  highest  frequency 
present,  namely  1.71  GHz.  All  the  methods  consistently  underestimated  the 
magnitudes  of  the  real  parts  of  the  last  two  pole  pairs.  The  Kalman  method  did 
not  help  with  this.  Apparently  the  highest  frequency  was  too  close  to  the 
Nyquist  frequency.  When  At  was  reduced  to  0.25  nsec  (Nyquist  frequency  * 

2.00  GHz),  the  results  improved  greatly,  and  the  errors  in  the  real  parts  of 
the  last  two  poles  were  random. 


APPLICATION  TO  EXPERIMENTAL  DATA 

Because  of  problems  with  the  tape  recorder  we  originally  planned  to  use, 
we  had  to  hand  digitize  the  experimental  curves  we  analyzed  and  punch  the 
values  on  cards  to  feed  into  the  computer.  The  oscilloscope  traces  were  en¬ 
larged  to  30.5  x  38.1  cm  (12  x  15  in.).  The  201  points  were  hand  digitized 
from  each  enlargement.  Each  set  of  data  was  interpolated  and  smoothed  to  512 
points  (at  equally  spaced  time  intervals) ,  using  a  cubic  spline  routine.  This 
procedure  was  performed  for  the  returns  from  the  hemisphere,  the  cylinder,  and 
no  object  (clutter). 

It  can  be  seen  from  Figure  8  that  the  returns  fircm  th*>  adjects  had  a  lot 
of  clutter  superimposed  on  them.  This  clutter  tat*??  $eres  waxfr  the  analysis 
of  the  return  to  find  the  poles.  If  the  cluttet  were  just  additive,  it  could 
be  removed  by  subtracting  the  object-free  return.  This  was  tried  but  was  not 
very  successful.  Any  reflections  arising  in  the  transmitting  circuitry  are 
convolved  with  the  impulse  response  and,  thus,  cannot  be  subtracted  out.  This 
apparently  was  the  problem  here.  It  might  be  possible  to  deconvolve  the  clut¬ 
ter,  if  enough  information  were  obtained  about  its  characteristics  to  allow 
distinguishing  it  in  some  way  from  the  signal. 

Tables  3  and  4  show  the  results  of  applying  the  various  methods  to  the 
cylinder  and  hemisphere  data.  The  total  data  span  for  each  was  10  nsec,  divid¬ 
ed  into  511  time  steps.  This  gives  a  data  spacing  of  0.01957  nsec.  The  theory 
only  applies  to  the  free-ringing  response  of  the  object,  so  the  first  point  to 
be  used  must  be  received  later  than  the  reflection  (if  any)  of  the  trailing 
edge  of  the  pulse  from  the  farthest  edge  of  the  object.  Since  the  objects  are 
each  12.7  cm  (5  in.)  in  diameter  and  the  pulse  width  is  about  1  nsec,  the  time 
from  the  beginning  of  the  reflection  to  the  first  point  to  be  used  is  about 
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i  +  2  x  12.7cm 
30cm/nsec 


■  1.85  mec 


Table  3,  Poles  from  Cylinder 

Nonlinear-Least  Squares 


Clutter  Not 
Removed 
(GHz) 

Clutter 

Removed 

(GHz) 

Clutter  Not 
Removed 
(GHz) 

Clutter 

Removed 

(GHz) 

-0.16 

+ 

j 

0.037 

-0.14 

±  j 

0.037 

-0.29 

+ 

j 

0.038 

-0.32 

±  j 

0.038 

-0.75 

t 

j 

0.22 

-0.24 

±  j 

0.11 

-0.54 

± 

j 

0.19 

-1.88 

±  j 

0.14 

-2.03 

± 

j 

0.38 

-0.82 

i  3 

0.21 

-1.24 

± 

j 

0.40 

-0.74 

4  j 

0.21 

-2.97 

± 

j 

0.65 

-1.09 

±  j 

0.31 

-3.09 

± 

j 

0.65 

-0.69 

±  j 

0.40 

-3.13 

± 

j 

0.82 

-0.68 

±  j 

0.37 

-2.33 

± 

j 

0.80 

-1.73 

±  j 

1.04 

-4.01 

± 

j 

1.02 

-6.11 

±  j 

0.82 

-2.16 

± 

j 

1.06 

-2.76 

± 

3 

1.40 

-1.56 

*  j 

1.56 

-3.45 

± 

j 

1.26 

-2.36 

± 

j 

1.41 

ible  3 

shows  the 

results  of 

applying  the 

linear  and 

nonlinear 

least 

squares  methods  to  the  cylinder  data  both  in  its  original  form  and  with  the 
clutter  subtracted.  Table  4  shows  the  results  from  applying  all  three  methods 
to  the  return  from  the  hemisphere  with  the  clutter  subtracted.  The  least- 
squares  methods  used  time  steps  of  0.03914  nsec  and  95  data  points.  The  Kalman 
method  used  time  steps  of  0.01957  nsec  and  200  data  points. 


Table  4.  Poles  from  Hemisphere,  Clutter  Subtracted 


Linear -Least  Nonlinear 

Squares _  Least  Squares _  Kalman 


-9.59  ± 

j 

0 • 1 1x1 0“18 

-0.29 

± 

j 

0. 

-0.60 

+ 

j 

0. 

-1.54  ± 

j 

0.17x10-14 

-8.33 

± 

j 

0.11x10“u 

-4.86 

± 

j 

1.23 

-1.41  t 

j 

1.12 

-0.80 

± 

j 

1.49 

-4.93 

± 

j 

3.32 

-1.72  ± 

j 

2.26 

-1.83 

± 

j 

2.19 

-1.28 

± 

j 

5.78 

-2.07  ± 

j 

3.57 

-0.07 

± 

j 

3.05 

-2.12 

± 

j 

8.08 

-1.81  ± 

j 

4.85 

-0.003 

± 

j 

5.03 

-0.92  ± 

j 

6.46 

-0.37 

t 

j 

6.37 

-2.55  ± 

j 

8.04 

-0.08 

± 

j 

7.66 

-3.94  ± 

j 

9.37 

-1.25 

± 

j 

8.56 
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It  can  be  seen  that  the  results  are  not  very  consistent.  Since  the  meth¬ 
ods  worked  well  on  synthetic  data  with  added  random  noise,  the  problem  here  is 
apparently  due  to  excessive  clutter  that  is  not  random  and  is  correlated  with 
the  signal.  Besides,  the  clutter  cannot  be  represented  as  a  sum  of  damped  ex¬ 
ponentials,  which  Prony's  algorithm  demands.  If  one  examines  the  clutter  care¬ 
fully  (Figure  8)  and  compares  it  with  the  impulse  generator  pulse  (Figure  5) ,  one 
sees  that  the  clutter  appears  to  be  a  series  of  reflections  or  miniature  repli¬ 
cas  of  the  impulse  generator  pulse.  We  tried  Prony's  algorithm  and  the  various 
methods  on  the  impulse  generator  pulse  as  received  by  the  TEM  horn  to  check  the 
tolerance  of  Prony's  algorithm  to  nonexponential  waveforms.  The  results  on  the 
impulse  generator  pulse  were  totally  meaningless.  Since  the  clutter  appears  to 
consist  of  miniature  replicas  of  the  impulse  generator  pulse,  we  would  expect 
that  Prony's  algorithm  will  not  work  well  on  such  a  waveform. 

To  obtain  meaningful  results  with  the  presently  available  techniques r 
the  clutter  must  be  reduced.  This  can  probably  be  achieved  by  using  fewer 
and  better  connectors.  Another  possible  way  to  deal  with  the  clutter  may  be 
to  deconvolve  it  by  homomorphic  deconvolution  techniques.11  This  may  work, 
but  only  if  the  frequencies  contained  in  the  clutter  are  vastly  different  from 
the  frequencies  in  the  returns. 

We  had  intended  to  do  runs  similar  to  those  described  above  on  the  return 
from  the  cone,  but  in  view  of  the  clutter  problems,  we  decided  to  forgo  the 
cone  runs  until  after  the  clutter  problems  are  solved. 


CONCLUSIONS 

The  time-domain  range  and  the  method  for  gathering  experimental  data  were 
described.  Prony's  algorithm  was  discussed  as  a  means  fot  processing  scattered 
returns  from  objects.  Improvements  in  the  algorithm  were  described  that  help 
it  to  better  process  noisy  synthetic  data.  It  did  not,  however,  do  well  on 
experimental  returns.  The  main  reason  was  clutter  that  cannot  be  represented 
as  a  sum  of  damped  sinusoidal  waveforms,  which  Prony's  algorithm  demands. 

The  needs  for  further  effort  are  as  follows: 

1.  Improve  experimental  techniques  so  that  the  data  obtained  are  in  a 
form  that  can  be  handled  by  presently  available  signal  processing  techniques. 
This  would  include  getting  rid  of  clutter  and  noise  as  much  as  possible.  Clut¬ 
ter  can  be  reduced  by  eliminating  unnecessary  connectors  and  replacing  those 
needed  with  connectors  that  are  good  at  high  frequencies  (15  GHz  and  up)  to 
reduce  reflections.  Also,  grounds  need  to  be  checked  to  eliminate  ground 
loops.  To  reduce  noise,  we  recommend  a  minicomputer  that  can  control  the  sweep 
of  the  sampling  scope.  In  this  way,  many  samples  of  the  same  point  can  be 
received  and  averaged. 
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2.  Improve  signal  processing  techniques  so  that  good  results  can  be  ob¬ 
tained  from  bad  (noisy,  cluttered)  data.  This  would  include  intensive  param¬ 
eter  studies  on  Prony's  algorithm  and  seeking  ways  to  deconvolve  clutter  and 
noise  from  the  sought  for  returns.  Also,  new  techniques  should  be  sought  and 
tried. 


REFERENCES 


1.  C.  E.  Baum,  On  the  Singularity  Expansion  Method  for  the  Solution  of 
Electromagnetic  Interaction  Problems ,  Interaction  Notes,  Note  88,  (11 
December  1971). 

2.  R.  Prony,  "Essai  Experimental  et  Analytique, "  J.  Ecole  Poly technique, 

(Paris,  1795),  2  (2),  pp.  24-76. 

3.  V.  C.  Martins,  J.  L.  Van  Meter,  J.  M.  Proud,  and  D.  J.  Fitzgerald, 
Picosecond  Pulse  Reflector  Antenna  Investigation,  RADC-TR-73-215 
(Griffiss  Air  Force  Base,  N.Y.,  July  1973). 

4.  A.  G.  Evans  and  R.  Fischl,  "Optimal  Least  Squares  Time-Domain  Synthesis 
of  Recursive  Digital  Filters,"  IEEE  Transactions  on  Audio  and  Electro¬ 
acoustics,  AV-21,  (February,  1973),  pp.  61-65. 

5.  J.  N.  Brittingham,  E.  K.  Miller,  and  J.  L.  Willows,  The  Derivation 

of  Simple  Poles  in  a  Transfer  Function  from  Real-Frequency  information, 
UCRL-52050,  Lawrence  Livermore  Laboratory,  (6  April  1976),  pp.  2-4. 

6.  R.  H.  Pennington,  Introductory  Computer  Methods  and  Numerical  Analysis, 
(New  York;  Macmillan,  1970) • 

7.  F.  B.  Hildebrand,  Introduction  to  Numerical  Analysis,  (New  York: 
McGraw-Hill,  1974) . 

8.  C.  L,  Lawson  and  R.  J.  Hanson,  Solving  Least  Squares  Problems 
(Englewood  Cliffs,  N.J.:  Prentice  Hall,  1974). 

9.  J.  L.  Junkins,  Optimal  Estimation  Theory,  NSWC  TO-K-45/74  (Dahlgren, 

Va.,  January  1975),  pp.  19-27. 

10.  M.  L.  Van  Blaricum  and  R.  Mittra,  ” Problems  and  Solutions  Associated 
with  Prony* s  Method  for  Processing  Transient  Data,*9  IEEE  Transactions  on 
Antennas  and  Propagation,  AP-26,  (January,  1978),  pp.  174-182. 

11.  S.  M.  Riad  and  N.  S.  Nahman,  "Application  of  the  Homomorphic  Deconvolution 
for  the  Separation  of  TDR  Signals  Occurring  in  Overlapping  Time  Windows," 
IEEE  Transactions  on  Instrumentation  and  Measurement,  IM-25,  (Decem¬ 
ber,  1976),  pp.  388-391. 


-rr 


\ 


tt: 


25 


DISTRIBUTION 


Naval  Air  Systems  Command 
AIR-31  OB 

ATTN:  J.  W.  Willis  (2) 

Washington,  DC  20361 

Commanding  Officer 

Naval  Ordnance  Station 

Gun  Systems  Engineering  Department 

ATTN:  Technical  Library  (Code  50D)  (2) 

Louisville,  KY  40214 

R.  C.  Hansen,  Inc. 

Box  215 

Tarzana,  CA  91356 

Naval  Postgraduate  School 
ATTN:  Michael  Morgan 

Monterey,  CA  93940 

The  Ohio  State  University 
Electro  Science  Laboratory 
1320  Kinnear  Road 

ATTN:  David  Moffatt  (4) 

Eric  Walton 
Columbus,  OH  43212 

Office  of  Naval  Research 
800  North  Quincy  Street 

Arlington,  VA  22217  (2) 

Commander 

Rome  Air  Development  Center 

ATTN:  Paul  Van  Etten 

Griffis  Air  Force  Base,  NY  13440 

Teledyne  Micronetics 
7155  Mission  Gorge  Road 
ATTN:  Steven  Weisbrod 

San  Diego,  CA  92120 

Director 

Naval  Research  Laboratory 
4555  Overlook  Avenue,  SW 
ATTN:  Philip  Moser 

Arthur  K.  Jordan 
Washington,  DC  20375 


DISTRIBUTION  (Continued) 


Office  of  Naval  Research 
Western  Regional  Office 
1030  East  Green  Street 

Pasadena,  CA  91106  (2) 

Electronics  Engineering  Department 
University  of  California 
Lawrence  Livermore  Laboratory 
Box  808 

ATTN:  E.  K.  Miller  (2) 

Livermore,  CA  94550 

Electromagnetics  Division 
National  Bureau  of  Standards 
ATTN:  Andy  Ondrejka  (4) 

Norris  S .  Nahman 
Robert  A*  Lawton 
Boulder,  CO  80303 

General  Dynamics 
Electronics  Division 
Box  81127 

ATTN:  Gus  Tricoles 

San  Diego,  CA  92138 

Systems  Applications 

Sperry  Corporate  Research  Center 

100  North  Road 

ATTN:  C.  Leonard  Bennett 

Sudbury,  MA  01776 

Effects  Technology,  Inc* 

Electromagnetics  Section 
5383  Hollister  Avenue 
ATTN:  Michael  L.  VanBlarium 

Santa  Barbara ,  CA  93111 

Picosecond  Pulse  Labs 
8663  Hollyhock  Lane 
ATTN:  James  R.  Andrews 
Lafayette,  CO  80026 

Department  of  Electrical  Engineering 
Texas  Tech  University 
ATTN:  Thomas  F.  Trost 

Lubbock ,  TX  79409 


DISTRIBUTION  (Continued) 


Department  of  Physics 
Catholic  University  of  America 
ATTN:  Herbert  Uberall  (2) 

Washington,  DC  20064 

Department  of  Electrical  Engineering 
University  of  Illinois  at  Chicago  Circle 
Box  4348  SEO-1104 
ATTN:  Wofrgang  M.  Boerner 

Chicago,  IL  60680 

Department  of  Electrical  Engineering 
Rochester  Institute  of  Technology 
One  Lomb  Memorial  Drive 

ATTN:  Tapan  K.  Sarkar  (2) 

George  Whitman  Reed 
Rochester,  NY  38677 

Department  of  Electrical  Engineering 
Michigan  State  University 
ATTN :  Kun-Mu  Chen 

East  Lansing,  MI  48824 

Department  of  Electrical  Engineering 
Virginia  Polytechnic  Institute 
ATTN:  Sedki  M.  Riad 

Blacksburg,  VA  24060 

Electromagnetic  Radiation  Analysis  Branch 
U.  S.  Environmental  Protection  Agency 
Office  of  Radiation  Programs 
Box  18416 

ATTN:  Richard  A.  Tell 

Las  Vegas,  NV  89114 

Electromagnetic  Sciences  Laboratory 
SRI  International 
ATTN:  David  M.  Bubenik 

Menlo  Park,  CA  94025 


DISTRIBUTION  (Continued) 


Defense  Technical  Information  Center 
Cameron  Station 
Alexandria ,  VA  22314 

Library  of  Congress 

ATTN:  Gift  and  Exchange  Division 

Washington,  DC  20540 

Local: 

E31  (GIDEP) 

E41  (Hall) 

F 

F10 

F12 

FI 2  (Hollmann) 

R 

X210 


(12) 

(4) 


(10) 

(6) 


